Draft Modified February 2, 2008 



Streaming Instabilities in Protoplanetary Disks 

Andrew N. Youdin & Jeremy Goodman 
Princeton University Observatory, Princeton, NJ 08544 

ABSTRACT 

Interpenetrating streams of solids and gas in a Keplerian disk produce a local, 
linear instability. The two components mutually interact via aerodynamic drag, 
which generates radial drift and triggers unstable modes. The secular instability 
does not require self-gravity, yet it generates growing particle density pertur- 
bations that could seed planetesimal formation. Growth rates are slower than 
dynamical, but faster than radial drift, timescales. Growth rates, like stream- 
ing velocities, are maximized for marginal coupling (stopping times comparable 
dynamical times). Fastest growth occurs when the solid to gas density ratio is 
order unity and feedback is strongest. Curiously, growth is strongly suppressed 
when the densities are too nearly equal. The relation between background drift 
and wave properties is explained by analogy with Howard's semicircle theorem. 
The three-dimensional, two-fluid equations describe a sixth order (in the complex 
frequency) dispersion relation. A terminal velocity approximation allows simpli- 
fication to an approximate cubic dispersion relation. To describe the simplest 
manifestation of this instability, we ignore complicating (but possibly relevant) 
factors like vertical stratification, dispersion of particle sizes, turbulence, and 
self-gravity. We consider applications to planetesimal formation and compare 
our work to other studies of particle-gas dynamics. 

Subject headings: hydrodynamics — instabilities — planetary systems: formation 
— planetary systems: protoplanetary disks 



1. Introduction 

A major uncertainty in theories of planet formation occurs embarrassingly early, during 
the formation of planetesimals. Collisions are unlikely to result in coagulation over a wide 
range of sizes, from mm to km, since available binding energies (chemical or gravitational) are 
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negligible comparable to kinetic energies (Chokshi et al. 1993; Youdin & Shu 2002, hereafter 
YS; Youdin 2003). Zero gravity experiments confirm destructive cratering during low velocity 
impacts (Colwell 2003). The hypothesis that planetesimals form by gravitational collapse of 
solids that settle to the disk midplane (Goldreich & Ward 1973; YS) remains controversial 
because it is uncertain whether protoplanetary disks are ever suitably quiescent or metal 
rich enough to allow gravitational instabilities to develop. 

The strong coupling between small solids and gas is both an obstacle to, and a necessary 
ingredient in, planetesimal formation theory. Forming planetesimals via gravitational insta- 
bilities would be trivial in a gas free disk, since collisional damping dominates viscous stirring 
in the absence of protoplanets (Goldreich et al. 2004a). Indeed Goldreich et al. (2004b) dis- 
cuss a second generation of planetesimals that could form in a gas free environment during 
the final stages of planet formation. However the first generation of planetesimals likely 
formed within the gas rich disks are observed to persist for several million years around 
low mass stars (Haisch et al. 2001). The outer planets of our solar system contain substan- 
tial amounts of this gas (Uranus and Neptune have several earth mass atmospheres, while 
Jupiter and Saturn are mostly gas) which presumably accreted onto solid cores assembled 
from planetesimals that, by this logic, must have formed in the presence of gas. While only 
trace amounts of gas are found on terrestrial planets, the simplest hypothesis is that inner- 
disk planetesimals also formed in the presence of gas. Indeed conditions for gravitational 
instability become favorable when the gas disk is only partly dissipated (Sekiya 1998, YS). 

The relative abundances of solids and gas are not yet tightly constrained by observations. 
At solar abundances the ratio of condensible solids to gas is of order 0.01, and is higher or 
lower depending on whether ices condense. The solid to gas surface density ratio, S p /S 9 , need 
not be fixed at solar abundances. YS and Youdin & Chiang (2004), hereafter YC, discuss 
and model enrichment mechanisms (e.g. radial migration and photoevaporation) that act to 
increase S p /E 9 . The ratio of space densities, p p /p g , is larger than S p /S 9 toward the midplane 
of a stratified disk. The extent of particle settling is limited by settling times (longer than disk 
lifetimes for sub-//m grains) and by turbulent diffusion, which can be generated by vertical 
stratification itself, among other possibilities. Assuming that Kelvin-Helmoltz instabilities 
trigger this stirring, Sekiya (1998) and YS showed that if S p /E 9 is augmented by a factor 
~ 10 the particle layer becomes so stratified that the gas and particle masses are equal in 
the particle sublayer, i.e. p p ~ p g throughout the layer. For even higher concentrations 
of solids, i.e. when the layer becomes particle-dominated, vertical shear instabilities are no 
longer capable of preventing particle settling and gravitational instabilities appear inevitable 
(YS). 

Given that interesting effects occur when the solid and gas densities become comparable, 
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Newton's Third Law tells us that we must consider the effects of frictional coupling on solids 
and gas equally. However most studies of midplane particle dynamics (see §7) do not fully 
treat the feedback of drag on the gas fluid. We consider a simplified model which treats the 
dynamics of both components self-consistently. We investigate the linear stability properties 
of a two fluid Keplerian disk, where a pressureless fluid represents particles of a specific 
size. It is well known that this system leads to steady state drift as angular momentum is 
transferred from the solids to the pressure-supported, and thus sub-Keplerian, gas (Nakagawa 
et al. 1986, see §2). The radial component of this drift globally redistributes solids on long 
timescales, leading to "particle pileups" since the accretion rate of solids is faster in outer 
disk (YS, YC). 

Here, however, we are concerned with local and more rapid consequences of orbital 
drift. By analogy with the two stream instability in plasma physics (Spitzer 1965), coupling 
between interpenetrating streams destabilizes linear waves. In our case the streams interact 
by drag forces, not electric fields. Our model does not include self-gravity. Nevertheless, 
unstable waves generate particle density perturbations. In principle these perturbations 
could be relevant to planetesimal formation, for instance by raising the particle density 
to a point where self-gravity induces collapse of the perturbations. We caution that the 
actual manifestation of particle-gas coupling may differ significantly from our model due 
to several simplifications. We ignore vertical structure in our background state, so our 
system is effectively an infinite cylinder and not a thin disk. Such an approximation may be 
warranted for vertical wavelengths smaller than disk scaleheights. Furthermore our model is 
linear, laminar, and inviscid, though a possible non-linear outcome of the instability is weak 
turbulence. 

This paper is organized as follows. Model equations and assumptions, for both steady 
state and perturbed motions, are presented in §2. Growth rates arising from the sixth order 
dispersion relation are numerically analyzed in §3. The relation between growth rates and 
wave speeds are studied in §3.3 by analogy with Howard's semicircle theorem. Eigenfunctions 
of a vertically standing waves are constructed in §4, allowing visualization of the fluid mo- 
tions. An approximate cubic dispersion that reproduces most features of the growing modes 
is derived in §5, allowing analytic investigation to complement the results of §3. Astrophys- 
ical applications considered in §6 include particle concentration (§6.1) and a comparison of 
growth rates to steady state drift (§6.2). We compare our work to other studies of dust layer 
dynamics in §7. A summary and conclusions are given in §8. 
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2. Basic Equations 

Our gas and particle "fluids" obey continuity and Euler equations for the evolution of 
particle (V p ) and incompressible 1 gas (V g ) velocity, here presented in a non-rotating frame: 



+ v-( Pp v p ) = 0, (i; 



dpp 
dt 

V-V 9 = 0, (2) 
<9V V — V 

^ + V P .VV P = -^ r -^_^, (3) 

^stop 

Ol Pg Cstop Pg 

where P is the gas pressure, p p and p g are the particle and gas spatial densities, respectively, 
and VLk oc r -3 / 2 is the Keplerian orbital frequency at cylindrical radius r. We ignore vertical 
stratification and self gravity for a simpler analysis, avoiding in particular the vertical settling 
and stirring of particles. The particle stopping time, t stop , is conveniently independent of p p , 
V p and ~V g , for the small particles (radius < 1 m at 1 AU) of interest prior to planetesimal 
formation. Epstein's law: 

Pg c 9 

applies when a < (4/9)A m f p , where a is the particle radius, A m f p is the gas mean free path 
(and A m f p ~ 1 cm at 1 AU), c g the gas sound speed, and p s denotes the material density of 
the solid. Particles larger than (4/9)A m f p , but small enough that the Reynolds number of 
the flow past the solid, Re = 4a|V p — V 9 |/(c 9 A m / p ) < 1 obey Stokes' law: 



stop 9o c A f ' 1 ' 

For generality we use the dimensionless stopping time parameter: 

T s = fi^-tstop (7) 

instead of referring to specific particle sizes and disk models. 

In this context, 2 fluid description of particle motions (as opposed to the kinetic theory 
approach) requires that solids be tightly coupled to gas. The criterion, t s « 1, ensures 



1 Since motions are very subsonic, this assumption filters sound waves from the analysis. 

2 A fluid description might also be possible given frequent intcrparticle collisions, but for small solids in a 
gas disk the stopping time is shorter than the collision time. 
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strong coupling to dynamical perturbations, while u;t s top <C 1, suffices for disturbances of 
arbitrary frequency, uj. We will not consider t s > 1. 

Since relative motions between solids and gas are slow compared to center of mass 
velocities (in equilibrium and for perturbations), we express equations (3, 4) in terms of 
relative motion, AV = V p — V 9 , and center of mass (henceforth COM) motion, V = 
(PpVp + p 9 V 9 )/p, with p = p P + p g the total density: 

<9V 

— + V • VV + F(AV 2 ) = -n 2 K v - VP/p , (8) 
^ + V • V(AV) + AV • VV + G(AV 2 ) = - A— + — . (9) 

(Jt Pg t sto p Pg 

The functions 

F(AV 2 ) = -V • (^AVAv) (10) 
G (AV 2 ) = ^AV • V (^AV )- ^AV • V f^AvV (11) 



can often be dropped due to the smallness of AV, in which case equations (8, 9) simplify 
considerably. (As the conditions for the neglect of these terms differs for equilibrium and 
perturbed motions, they will be addressed subsequently.) Another advantage of this formu- 
lation is that drag forces do not appear in COM evolution, and gravity, including self gravity 
were it included, is absent from the relative motion equation. 



2.1. Steady Drift Solutions 

The time steady, axisymmetric drift motions of solids and gas have been well studied 
(Nakagawa et al. 1986), and are simple to rederive from equations (8, 9) with F = G = 0: 



U = W = AW = (12) 



2 , r i)P ( , Pg. 
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where V = Ur + V(f)+Wz, AV = AUr+AVfi+AWz, overbars denote steady state solutions, 
Vk = &K r is the Keplerian circular speed, and 

1 dP ( c \ 2 

^-2^qm^~\t K ) (16) 

measures the radial pressure support. In standard models of planet forming disks, rj ~ 10~ 3 
at 1 AU. Neglect of the non-linear drift terms (F, G) is, in practice, always justified for 
equilibrium solutions. These terms would only merit inclusion in the unlikely scenario that 
pressure strongly dominated gravity, r] >> 1. 

Our equilibrium solutions (12 - - 15) contain no vertical motion, due to the neglect 
of vertical gravity. The center of mass is fixed in radius, orbiting as a gas supported by a 
pressure P, but with a mean molucular weight augmented by pj p g . For outwardly decreasing 
pressure, i] > 0, sub-Keplerian gas robs particles of angular momentum, giving inward 
migration of solids, U p = (p g /p)AU < 0, and outward migration of gas, U g = — (p p /p)AU. 
As it will set the scale for wave speeds, we introduce the unweighted sum of the radial drift 
speeds: 

U =u I U - 2p p9 ~ Pp vVkTs (17) 

u sum -u p + u g - z Pg p2 1 + {TsPs/p)2 {U) 

which has an interesting density dependence. In the test particle limit, p p — > 0, U sum — > 
Up — > —2t]Vkt s . With increasing particle concentration, the gas migrates faster at the 
expense of the solids, until U sum = for equal densities. At p v = 3p g , U sum = t]Vkt s /A 
reaches its (positive) maximum. For even larger particle concentrations U sum declines as the 
gas pressure is diluted by the particle mass. 



Azimuthal drift is much slower than radial, AV / AU = p g r s / (2p) for tight coupling. For 
loose coupling, AV — > tjVk as particle and gas trajectories approach unperturbed Keplerian 
and pressure supported rotation, respectively. The previous caution about applying fluid 
equations for r s ^> 1 can be ignored for these equilibrium solutions provided the disk varies 
negligibly over a stopping length t stop AU ~ r/r r. Finally note that our unstratified 
model, with dV/dz = 0, ignores the vertical shear, which can generate Kelvin- Helmholtz 
instabilities. 



2.2. Localized Perturbation Equations 

We consider perturbations to steady state motion (equations 13 — 15) on length scales 
much shorter than disk radii (indeed, shorter than the thickness of the gas layer, H g ~ 7] l ^ 2 r). 
This allows a local treatment, with Cartesian coordinates corotating about a fixed radius, 
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r a , with the orbital frequency of the local COM, f2 D = V t f,{r )/r . In the new coordinates: 

x = r — r Q , (18) 
y = r o (0 - ft t) , (19) 

we approximate differential rotation, as usual, by plane parallel flow with linear shear, V = 
—qQ D xy, where q ~ 3/2 for the nearly Keplerian profile. Within this approximation drift 
motions are radially constant, dAV/dx = 0. We decompose fluid variables into steady 
backgrounds and perturbations: 

V = -lQ xy + v(x,z,t), (20) 

AV = AUx + ~A~Vy + Av(x, z, t) , (21) 

p = Po [l + 5(x,z,t)}, (22) 

P = Po[-g e x + h(x,z,t)], (23) 

where g e = —dP Q / dr\ To / p Q = 2rjVl 2 r p g / p Q > 0. Henceforth we drop overbars and the 
subscripted o's from unperturbed states. Perturbations are axisymmetric, which avoids 
stretching by radial shear, and are given a Fourier dependence: 

f(x, z, t) = fexp[i(k x x + k z z - ut)] , (24) 

with real wavenumbers, k x and k z , and a complex frequency, 

uj = + is , (25) 

with a wave frequency, u^, and growth (or damping) rate, s. 3 

The linear perturbation equations read: 

— uvv — 2Vtvx + —uy + F' = — ikh — 5g e x (26) 

^* - . ., Arr ^/ (AV + 5AV) , /l 

-iujAv - 2£lAvx H — Aim/ + ik x AUv + G = -- — + ik— (27) 

2 Jgtstop Jg 

-iu5 + ik-v = (28) 
kv = f p k-Av + f g k x AU5 (29) 

and the terms that arise from perturbations of equations (10, 11) are: 

F' = if g [(f p k- Av + f g k x AU5)AV + f p k x AUAv] (30) 
G' = ik x AU[(f g -f p )Av-f g AV5] (31) 



3 Modes with s > and uj^ ^ arc often called "overstable" to distinguish them from non-oscillatory 
instabilities. 
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where v = ux + vy + wz, Av = Aux + Avy + Awz, f g = p g /p is the equilibrium gas fraction, 
and the particle fraction, f p = 1 — f g . 

Equations (26 — 31) define a 6th order (one less than the number of time derivatives 
because of the incompressibility constraint) dispersion relation for u, whose solutions we 
investigate in §3. Equivalent results are obtainted by perturbing particle and gas equations 
(1 — 4) directly, but the relative velocity formulation allows analytic simplifications. For low- 
frequency waves with u < Q, it is safe to neglect F' and G', but all terms are included in our 
numerical solutions. In §5, further approximations allow us to derive an approximate cubic 
dispersion relation that describes the growing modes while filtering three strongly damped 
modes. 



3. Growth Rates and Wave Speeds 

We investigate the growth rates and wave speeds, f wa vc = ^n/k x of linear disturbances. 
The eigenvalue problem for u/Q (prescribed by equations 26 — 29) is uniquely specified 
by four quantities: the stopping time parameter, r s = Qt stop , and the equilibrium solid-gas 
density ratio, p p /p g , which define the background state, and two normalized wavenumbers, 
K x = r)rk x and K z = rjrk z . Of the six modes, three decay within a stopping time, and are of 
little physical interest. The other three modes, of which two are modified epicycles and the 
other is a uniquely two fluid secular mode, can be slowly damped or growing. The secular 
mode gives the fastest growth. This section investigates the fastest growing waves, consid- 
ers whether turbulent viscosity can stabilize short wavelength modes, and relates growth 
rates and wave speeds to the background flow, by analogy with Howard's (1961) semicircle 
theorem. 



3.1. General Features 

Figure 1 (left) plots growth rate contours versus stopping time and p p /p g - We hold 
K z — 1 fixed (more on this choice later) and maximize the growth rate with respect to 
K x , which is plotted at right. Growth is possible for all values of r s and p p /p g , on slower 
than dynamical timescales. In the well coupled regime, t s <C 1, peak growth rates increase 
as s oc t s (for fixed p p / p g ) since looser coupling increases the relative motion needed for 
instability. The growth rates peak for marginal coupling and decrease for r s > 1, but we 
ignore this regime where a fluid description of the solids is questionable. 

The density dependence in figure 1 is more complicated, and has nothing to do with self- 
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Fig. 1. — (Left) Contours of growth rate (times orbital period) vs. stopping time and solid to 
gas density ratio. The vertical wavenumber is fixed at K z = 1 while K x is varied to maximize 
growth. (Right) Radial wavenumbers (solid contours) of these fastest growing modes. The 
power law fits (dotted contours) follow equation (32). See text (§3) for discussion. 



gravity, which is not included. As expected, growth rates decrease in the test-particle limit, 
as p P l p g — > 0, and as p v j p g — > oo when solids are unaffected by drag. More surprisingly, 
growth is diminished in a narrow region around p p = p g . We will show that this is related to 
vanishing wave speeds. Consequently, two lobes of relatively fast growth (around p v j p g fa 0.2 
and 3) exist where particle gas feedback is significant, but not so close to p p = p g that waves 
stagnate. The fastest growth occurs in the particle-dominated lobe. 

The K x values of figure 1 (right) are well fit by a broken power law: 

K _ f Vtj;)-^ if f g < 1/2 

X ~ { V2r-^f-^ if /, > 1/2 [d2) 

except for the spike at f g — 1/2. The preferred radial wavelengths decrease with the particle 
stopping time as K x oc r s 1 ^ 2 . The increase in K x with p p / p g is related to the density 
dependence of the gas stopping time, t st0 pPg/Pp (see equation 4), and stopping time for 
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Fig. 2. — {Left) Contours of growth rate, log 10 (s/Jl) (solid lines), and damping rate, 
log 10 (— s/Q) (dashed lines), vs. K x and K z for p v j ' p g = 0.2 and t s = 0.01. Two regions 
contain the fastest growing modes: K z > 10 2 , K x tv 80 (the darkly shaded region) and along 
K z ps r s Klf^ (the dotted line in both figures). (Right) Wave speed, f wav e = wm/k x , in units 
of —i]Vkt s , for the same modes. The contours from 0.1 to 1.1 increment by 0.2 and s peaks 
along this gradient in phase speed. Darkly shaded, large phase speed regions to the right 
(and upper left) correspond to damped (and very weakly growing) modes, respectively. 



relative motions, t stop p g /p (see equation 9). The two fluids become more tightly coupled as 
particle concentration increases, resulting in shorter wavelength growing modes. 



3.2. The Long and Short of It 

Figure 2 (left) plots growth rate versus wavenumbers (K X ,K Z ). Specific values of the 
density ratio and stopping time (p p /p g = 0.2, t s = 0.01) are chosen, but the qualitative 
features are rather general. We can identify two "ridges" along which growth rates peak. 
The short wavelength branch follows vertical contours (K x ps 80 and K z > 100 in the 
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figure) while the long wavelength ridge falls diagonally along K z ~ t s K% (a generalization of 
equation 32 for K z 7^ 1 that ignores the density dependence). A smooth transition between 
the two ridges occurs around K x ~ K z ~ l/r s . Very long wavelength modes (K x < 1, 
K z < r s ) are damped by frictional dissipation and angular momentum gradients. Much of 
this paper chooses the long wavelength branch by fixing K z — 1. The exact value chosen 
is arbitrary as growth rates vary little along this ridge. This subsection analyzes the short 
waves which have larger growth rates, but could be less relevant if turbulent viscosity were 
included in the analysis. 

3.2.1. Short Wavelength Limit: K z > K x 

The short wavelength behavior is described by a dispersion relation that is independent 
of K z for K z ^> K x , as seen in figure 2 . Physically, radial pressure perturbations become 
negligible. Figure 3 plots the maximum growth rates, and the K x values of the fastest 
growing modes, in this large K z limit. For a gas-dominated system, the short waves grow 
marginally faster than the long wave branch, the difference is less than a factor of 5 in the 
gas dominated regime, p v j p g < 0.2. The preferred radial wavenumber is nearly constant, 
K x ~ 1/t s , in the gas-dominated regime. The growth rate skyrockets as the density ratio, 
Pp/Pg, increases toward and above unity. This contrasts with the behavior we saw in the 
long wavelength case where growth is suppressed near equal densities. The K x value that 
maximizes growth increases with particle fraction, and does not keep the characteristic value 
l/r s . 

The real frequencies of these modes (in both the particle- and gas-dominated regimes, 
not plotted) are near, but slightly below, the dynamical frequency. Longer wavelength modes 
have lower frequencies. Indeed, figure 2 shows that f wav c = oj^/k x is similar for the fastest 
modes at all wavelengths. Hence modes with higher k x must have higher frequencies. The 
terms in equation (30, 31) must be kept to describe the short wavelength, high frequency 
modes. 



3.2.2. Turbulent diffusion 

Another factor to consider in determining the relative significance of short or long waves 
is turbulent diffusion. Turbulence is not included in our model, and so we cannot be certain of 
its effects. If it introduces viscous diffusion, short wavelength modes would be preferentially 
damped. To estimate the relevance of this effect, consider the diffusive timescale, t D ~ 
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Fig. 3. — Maximum growth rates (solid line) and fastest growing radial wavenumber (dashed 
line) versus solid to gas density ratio for r s = 0.01 in the limit K z ^> K x . The growth 
rates are below the upper limit implied by the semicircle theorem (dotted lines), except for 
a narrow region near equal densities. 



Considerable uncertainly surrounds the appropriate value for a. The values invoked to 
explain accretion onto young stars, ranging from 10~ 4 — 10~ 2 at least, may not apply here. 
Even if accretion is driven by turbulent diffusion, disks likely contain spatial inhomegeneities 
(disk midplanes which may be more quiescent) and experience temporal evolution (accretion 
rates decrease with age). 

It is more relevant to consider diffusivities needed to maintain a thin, but finite density, 
particle layer. A simple balance between settling and diffusion (see e.g. YC) for well-coupled 
particles (r s <C 1) suggests that a ~ r]T s [H p / (r/r)} 2 is required for sedimentation to a particle 
scaleheight H p (which is assumed to be thinner than the gas scale height). If H p m rjr (thinner 



47r 2 /(k 2 D) with the diffusion coefficient parametrized by the usual prescription D = ac 2 /Q. 
Growth outpaces diffusion if sip > 1, or equivalently: 




(33) 
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layers may be strongly Kelvin-Helmoltz unstable), equation (33) gives K < 2ir^ ' s/ \VLt s ). 
Since s < Qr s (for p p < p g at least), short modes with K 1 should be strongly damped. 
Even longer wavelength modes, with K ~ 1 — 10, are affected viscous diffusion, according to 
this analysis. However, viscous effects can sometimes destabilize disks (Schmit & Tscharnuter 
1995). This issue merits further study. 




Fig. 4. — {Left) Growth rate, \og w {s/Q), {solid contours) and decay rates, log 10 (— s/Q) 
{dotted contours) vs. solid to gas density ratio and radial wavenumber for t s = 0.01 and 
K z = 1. Two lobes of rapid growth are centered on p p /p g ~ 0.2 and 3, with suppressed 
growth near p p = p g . {Right) Radial wave speed contours in units of t]VkT s for the same 
modes. The phase speed changes sign across p p = p g . The dashed contours in both plots 
indicate the location where v wave = U snm /2. 



3.3. Phase Speeds and the Semicircle Theorem 



Figure 2 {bottom) plots the wave speed, f wav e, of the modes whose growth rates were 
shown in figure 2. For the mode of interest, we see two plateaus of nearly constant phase 
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speed. 4 The steep transition between these values overlaps the ridge of large growth rates 
in figure 2. 

These results are analogous to Howard's semicircle theorem for the Kelvin-Helmholtz 
instability, which we summarize briefly (see Kundu & Cohen (2002) for a derivation). Howard 
(1961) found that one dimensional modes, oc exp(st — + ik x x), have a wave speed that 
lies between the minimum and maximum speeds in the shearing flow, V m \ n < uj^/k < V max . 
The semicircle theorem: 

[ux/k - l(V max + V min ) 2 ] 2 + {s/kf < [l(V max - V min ) 2 ] 2 (34) 

says that the complex wave velocity of an unstable mode lies in a semicircle (since only 
positive growth rates are considered) of radius V max — V min . This imposes a limit on the 
maximum growth rate: 

s — 2(^ / i nax ^min)j (35) 

which is achieved for a phase speed midway between the allowed range. 

The physical differences between our streaming instability and the Kelvin-Helmholtz 
instability cannot be overstated: two interpenetrating, unstratified, rotating fluids with 3D 
velocities and 2D wavenumbers versus a single, plane-parallel, neutrally buoyant fluid with 
2D velocities and ID wavenumbers. It is remarkable then that our modes behave as if a 
modified semicircle theorem applied to them. As no proof of an analogous theorem for our 
problem exists, we describe the similarities. First, the radial phase speeds of our secular 
growing mode fall in the range, < |f wavc | < |t/ SU m|, 5 where U sum is the sum of particle 
and gas drift velocities, see equation (17). As U sum is positive (negative) in a particle (gas) 
dominated layer, respectively f wa ve has the same sign as U sum . Secondly, growth rates are 
largest U sum /2 and vanish near the endpoints of the allowed range. Since 

U sum = for equal densities (p p = p g ), the semicircle theorem is consistent with the finding 
that growth is weakened in this case (though not in the large wavenumber limit, as we have 
discussed). The peak growth rates are indeed bounded by s < \k x lI sum \/2 except for a very 
narrow region around p p = p g . 

To demonstrate the generality of these findings, figure 4 plots the growth rates and wave 
speeds, versus K x and p p /p g - At small K x , wave speeds, i> wave , approach U sum , for instance 
with p p / p g = 0.1, f W avc — ► U sum ~ —1.5r)VKT s as the contour indicates. For large K x , 



We can ignore the discontinuities in the upper left corner and far right hand side of the plot, which 
correspond to different, epicyclic, roots that happen to give larger growth rates in this region of phase space, 
which is generally uninteresting as growth rates are small. 

5 For p p > p g , the behavior is a bit more complicated. This case will be discussed shortly. 
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^wave — > (ignoring the mode switching of the dark region in the upper left corner). Growth 
rates are largest roughly midway through this transition, where f wav e ~ U smn /2, as indicated 
by the dashed contours. The suppression of growth for p p = p g , when t> wavc ~ U sum = and 
the "radius" of the semicircle vanishes, is clear. 

The transition in v wscve from U sum to with increasing K x has an added wrinkle in the 
particle-dominated case. The wave speed first rises slightly above U sum (which is clearly not 
a strict upper limit) before dropping to zero, as can be seen by following the contours in 
figure (4, right) for p p /p g > 1. The analogy to the semicircle theorem is still relevant, as the 
fastest growth occurs for f wav c ~ U snm /2. 

The analogy to the semicircle theorem connects the background flow to wave proper- 
ties. The free energy of interpenetrating streams undoubtedly plays a role, but from this 
perspective, it is surprising that the velocity scale is set by the sum rather the difference of 
radial streaming velocities. More study of two coupled, rotating fluids should increase our 
understanding of this system. 



4. Eigenfunctions: Fluid Motions and Density Perturbations 

Having investigated the eigenvalues, i.e. the growth rates and phase speeds of a mode, we 
consider the eigenfunctions, i.e. the Fourier amplitudes, /, that give us the fluid variables 
via equation (24). An individual mode has a vertical phase speed, u)&/k z , which can be 
eliminated by linearly superposing pairs of modes with opposite signs of k z . Under a vertical 
parity transformation, the vertical velocity is odd, while u and all other Fourier amplitudes 
are even. The vertical standing waves have forms: 

/odd = U(if exp[i(k x x - ut)}) sin(k z z) (36) 

/even = $l(f exp[l(k x X - Ut)]) COs(k z z) (37) 

for the odd (vertical velocity) and even (all other) variables, respectively. 

Figure 5 plots particle velocities for a rapidly growing mode. The K x value gives the 
fastest growth rate for K z = 1 as in figure 1. Since vertical wavelengths are longer than 
radial, tight coupling of particles to the incompressible gas, du g /dx + dw g /dz = 0, causes 
vertical velocities to dominate. Recall from figure 1 (right) that elongation of the fastest 
growing modes is more (less) pronounced for tighter (looser) coupling. The vertical velocities 
flow toward density maxima for this growing secular mode. The pair of epicycles are weakly 
damped (s m —2.4 x lCT 3 f2 and s m —8.9 x lCT 3 f2) for these parameters. Their flow 
patterns are similar to figure 5 except vertical velocities flow toward density minima. Gas 
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Fig. 5. — Instantaneous (perturbed) particle velocity, v p , in the x — z plane with a greyscale 
image of azimuthal velocities (white is positive) for a growing mode with K x = 5, \K Z \ = 1, 
t s ps .044, p p / p g = 0.2. Gas velocities are very similar due to strong coupling. The density 
is very nearly in phase with the azimuthal speed, so the vertical flow is channeled to high 
density regions. The ratio of azimuthal to vertical velocity amplitudes is |v p |/|w p | ~ 0.66. 
The radial to vertical ratio, |-u p |/|w p | ~ K z /K x = 0.2, follows from near incompressibility. 
This mode has a growth rate s/Q ~ 2.9 x 10~ 3 and a phase speed, oJ^/k x = — 0A2\AU\. 



and particle velocities are not well coupled for the three strongly damped modes. The gas 
is nearly stationary so particle motion leads to rapid decay. 

Figure 6 shows perturbed relative velocities for the same mode as figure 5 with the 
perturbed density in greyscale. This relative motion is predominantly radial, even accounting 
for azimuthal velocities. The correlation of density with Au can be derived from continuity 
equation using s. 

These eigenfunctions cannot be fit into a finite thickness dust layer. This is clear from 
equations (36 and 37) and figure 5 which show that vertical velocities are maximized where 
density (and other component of velocity) vanish, and vice versa. A more complicated model 
that includes either stratification or a free surface between the particle layer and overlying 
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Fig. 6. — Perturbed relative motion of solids and gas, Av, for the same mode as figure 
5. The greyscale image shows density perturbations (white is positive). The radial relative 
motion dominates the azimuthal, |Aw|/|Am| 0.15, and vertical, |Aw|/|Am| 0.11, speeds. 
Density perturbations correlate with relative motion. 



gas layers, would give eigenfunctions with more realistic boundary conditions. 



5. Terminal Velocity Approximation 

A simpler description of unstable modes, which filters the three strongly damped roots, is 
achieved by assuming that relative velocities reach "terminal velocity," so that drag forces ad- 
just quasistatically to pressure forces. This amounts to neglecting all terms on the left-hand 
side of equation (9), both in equilibrium and in perturbation. Thus AV = — (VP/ p)t stop , 
and perturbations obey 

Av + 5AV ^ ikht stop . (38) 

This approximation, which ignores inertial accelerations, is valid so long as K <C l/r s and 
r s < 1. 
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A cubic dispersion relation: 

Q S --Q 2 (^ +2 ^) - © (§)W*(§)V />.=<> (») 

results from equations (26, 28, 29, and 38), see the appendix for intermediate equations. 
The roots of this cubic reproduce the results of the full system of equations to very good 
approximation when the stated assumptions are met. 



5.1. Stability of Inplane Motions 

When K z = 0, equation (39) gives: 

-iu/SI = -f p r s - 2if g K x r s (40) 

so that all modes are damped, s < 0. The full equations also lack growing modes for k z = 0. 
Fluid motions in this case are quite limited, especially given gas incompressibility. Equations 
(26, 27) show that w = Aw = 0. 6 The gas incompressibility condition, k x u g = 0, requires 
u g = for a non-trivial mode. Thus gas velocities are one dimensional, azimuthal. We must 
allow two dimensional waves, and three dimensional motion, to find secular instability. 



5.2. Equal Mass 



When the mass densities of solids and gas are equal, so that f p — f g — 1/2, the constant 
term in (39) vanishes. In this case we have a static mode, uj = 0, and the quadratic roots: 
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(41) 



The growth rate s < for for any (real) choice of parameters. 7 Notice that as t s — > the 
modes approach the modified epicyclic frequency, uj/Q — > ±K Z /K (this is the frequency of 
inertial oscillations in a single incompressible fluid). Numerical solutions of the full equations 
give suppressed, but actually non-zero, growth when densities are nearly identical. Also 
small- wavelengths (K > l/r s ) behavior, in which growth rates increase near equal densities, 
see figure 3, is not captured in the terminal velocity approximation, as terms giving small 
scale accelerations have been dropped. 



6 If s = —l/(fgt s top), then Aw ^ is possible, but this damped mode is not of particular interest. 
7 This follows trivially from the fact that 3?(y / (1 + ia) 2 — b 2 ) < 1 for all real values of a and b. 
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5.3. 



Series Solutions 



Aside from the special cases above, it is more enlightening to consider solutions to the 
approximate dispersion relation, equation (39), as a series expansion in r s <C 1: 



(Series solutions of the full sixth order dispersion relation have been done, but are not 
presented here.) Two of the three modes described by equation (39) are epicycles (inertial 
oscillations) to leading order, with ujq/VL = ±K Z /K. The first order correction, ui/Q — 
—f g fpK x — if p (K x / K) 2 /2 shows that these modes are damped to lowest order. With the 
full set of equations, epicycles can grow for K z ^> K x , but growth rates of the secular mode 
are always faster. 

The third, secular root is oscillatory to leading order: 



Thus the leading order wave speed is uj\T s /k x = U sum , the sum of the equilibrium drift speeds 
of gas and solids, see equation (17). This agrees with the wave speed of growing modes for 
small K x , before higher order corrections become significant. Since 0J2 = for this mode, 003 
gives the leading order growth rate: 



This rate is always positive, but nominally third order in t s <C 1. However growth rates are 
larger for K x » K z . The growth is maximized at K z t s K 2 as a higher order expansion 
(and e.g. figure 2, top) shows. This asymmetry in wavenumbers makes the growth rate first 
order in t s . 

These low order expansions confirm some basic results about growth rates and wave 
speeds. Most importantly, we demonstrate that the secular mode is responsible for fastest 
growth, not the pair of epicycles. 



LU = LUq + UOiT s + U0 2 T 2 + ... . 



(42) 




(43) 




(44) 



6. Applications 



6.1. 



Particle Concentration 



Two fluid instabilities like ours could aid planetesimal formation by generating particle 
density perturbations. With the aid of self-gravity these perturbations could eventually 
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Fig. 7. — Contours of A$/A h , the ratio of particle density perturbations to radial pressure 
gradient perturbations, in r s — p p /p g space for the growing modes of figures 1. Density 
perturbations are sizeable for p p / p g < 1. The largest relative density perturbations, around 
Pp/pg ~ 1, correspond to slowly growing modes. 



collapse to solid densities. Perturbation amplitudes are arbitrary in a linear analysis, so 
inferences about non-linear development are speculative. To estimate the relevance of density 
perturbations, we compare the perturbation amplitudes of particle density, As = \S\/f p and 
radial pressure gradients, Ah = \k x h/g e \. Figure 7 shows that As > Ah, suggesting that 
density perturbations are significant. By comparison with figure 1, we see that the regions 
of largest growth rates do not coincide, and are somewhat anti-correlated actually, with the 
largest density perturbations. 

We briefly justify using radial pressure gradients as the scale to compare the density 
perturbations. Vertical gradients of pressure perturbations are smaller (since k z < k x ) and 
more importantly have no background value in our unstratified model. Perturbation ve- 
locities have six components and can be compared to several different background speeds, 
including drift, Keplerian shear, and tjvk, the amplitude of pressure supported sub-Keplerian 
rotation. However the amplitudes \v\/t)Vk and |Av|/A£7 are similar to Ah, i.e. somewhat 



smaller than Ag. 

The mass of solids in an unstable mode varies considerably over the wide range of 
possible wavenumbers. Let us conservatively take a small-scale mode with K x ~ K z ~ 100, 
in which the solid mass is 

£ p 27r27r27r 10 20 
M k = -JL — — — ~ g, (45) 



Hp k x k z ky 

where the particle surface density S p ~ 10g/cm 3 , 



kyTjr 

i] ~ 10~ 3 , and the particle scale height, 



H p rjr, is the value for stirring by Kelvin Helmholtz instabilities (refs). This is thin enough 
so that p p > O.lpg. Since our modes are initially axisymmetric, k y r]r < 1 is possible, but even 
if azimuthal breakup occurs on scales comparable to the radial wavelength, k y r\r ~ K x ~ 100, 
M k is comparable in mass to a 100 km planetesimal. Thus while non-linear development 
is unclear, the density perturbations induced by streaming instabilities contain more than 
enough mass to make healthy planetesimals. 
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Fig. 8. — Growth rates and equilibrium drift rates of solids and gas versus the density ratio, 
pp/ p g for f] = 2 x 10~ 3 ,t s = 0.01. The growth rates are faster than the drift rates for the 
given values. This is true more generally as well, see text. 
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6.2. Growth vs. Drift Rates 

Our local treatment of the instability is valid only if the growth is faster than global 
disk evolution, including changes in surface density or temperature. Single fluid accretion 
disk models evolve on timescales > 10 5 — 10 6 years, while passive disks evolve more slowly. 
Global redistribution of solids, at the equilibrium radial drift speed U p , changes particle 
surface densities (YS, YC). Gas densities are less subject to change because drift speeds are 
smaller (when p p < p g ) and because drift rates are much smaller for the majority of the gas 
mass, which lies outside the particle-dense midplane. 

To justify the local treatment of the instability, we compare growth rates to equilibrium 
radial drift, which leads to a global redistribution of solids (YS, YC). Figure 8 shows that 
growth rates are at least an order of magnitude faster than the particle drift rate, U p /r. Gas 
drift rates, also slower than growth rates, are shown as well. In a stratified disk, the gas drift 
rate will decrease with height as the particle concentration drops, while the particle drift 
rate asymptotes to a constant value. For all reasonable values of stopping time and pressure 
support, f], growth should still dominate. Both the growth and drift rates are linearly 
proportional to stopping time for t s <C 1. A hotter disk, i.e. larger rj ~ c 2 g /v 2 K oc T, has a 
faster drift rate, while growth rates are unaffected by rj. 8 As long as disks are not hotter 
than standard models by more than an order of magnitude (a safe assumption according 
to observations), we conclude that growth rates are robustly faster than drift rates. This 
enhances our confidence that the instability can be astrophysically relevant. 



7. Other Work 

There have been several previous studies of the linear stability of gas-and-dust mixtures 
in the context of the protosolar nebula. These works include physical effects that we have 
neglected, such as self-gravity and vertical stratification. But because of various restrictions 
on the types of perturbations considered, none found the modes described in this paper. 
Coradini et al. (1981) and Noh et al. (1991) wrote down two-fluid equations including self- 
gravity but permitted only horizontal motions. Sekiya (1983) allowed vertical motions but 
worked in the tightly-coupled limit where the velocity difference between the two fluids is 
negligible. All of these authors found gravitational instabilities at sufficiently high densities. 
Ishitsu & Sekiya (2003) and Garaud & Lin (2004) examined the stability of the vertical shear 
between the settling dust layer and the overlying gas but also neglected the slippage between 



The fastest growing wavelength is shorter in a hot disk, with k oc rj. 
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the two fluids. 

Like the present paper, Goodman & Pindor (2000) found an instability driven by drag 
rather than self-gravity, but theirs was not a complete two-fluid analysis. Following Goldre- 
ich & Ward (1973), GP treated the dust as a monolithic though dilute layer, the drag being 
exerted at its top and bottom surfaces by boundary-layer turbulence driven by the difference 
in orbital velocity between the dust-laden and dust-free components. The dynamics of the 
gas was not treated explicitly, its effects on the dust layer being parametrized in terms of 
the orbital velocity difference. GP did however emphasize that they expected the existence 
(though not the growth rate) of drag instabilities to be independent of many physical de- 
tails provided that the drag is a collective effect: in other words, that the drag on a dust 
particle depends upon neighboring particles and is not linearly proportional to dust mass. 
In GP's case, the collective property derives from the assumption that the drag depends 
on conditions at the surface of the dust layer only, so that the average drag per particle is 
inversely proportional to the column density of the dust layer. We, however, have resolved 
the vertical dimension explicitly, so that the drag on a small dust particle depends on its 
motion relative to the local gas only. But because of the backreaction of the dust on the gas, 
the drag is collective: as the local ratio p p /p g increases, the relative velocity and the drag 
per particle decrease. Apparently, this is enough to support an instability despite the many 
simplifications assumed for our background state. 

8. Conclusion 

We describe a two-fluid streaming instability relevant to protoplanetary disks of particles 
and gas. Unstable modes are powered by the relative drift between the two components, a 
universal consequence of radial pressure gradients. Rotation, which is Keplerian in our model, 
is another necessary ingredient. Growth occurs despite the fact that the two components 
interact only via dissipative drag forces. The robust instability has growth rates slower than 
dynamical times but faster than drift times. The fluid motions generate particle density 
enhancements, even in the absence of self-gravity, which could trigger planetesimal formation. 
The physics of our disk model was simplified considerably (see §2). Numerical studies that 
include vertical stratification, a dispersion of particle sizes, and non-linear effects could 
elucidate the role of such instabilities in protoplanetary disk evolution. 

We are grateful to Andrew Cumming, Marc Kuchner, Doug Lin, Gordon Ogilvie, and 
Scott Tremaine for helpful comments, many of which were made during the successful KITP 
workshop on Planet Formation. This material is based upon work supported by the National 
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of space science. This research was supported in part by the National Science Foundation 
under Grant No. PHY99-07949. 



A. Simplified Equations of Motion 

Derivation of the cubic dispersion relation, equation (39), uses the terminal velocity 
approximation, equation (38), with equations (26, 28, and 29) and neglecting the F' and G' 
terms. Directly solving for the characteristic equation of this set introduces terms, including 
a quartic in u, which must be dropped for consistency with the r s C 1 approximation. 
Alternatively, elimination of pressure and relative motion variables gives the following three 
equation set: 

(cu-f g k x AU)5 = 2f p k x r s v (Al) 
-iu;r-2Qk z v = k z ttAU5/r s (A2) 

(-^ + § fpTsQ ) V + W^ r = -W f ° QAm (A3) 

where T = k z u—k x w, is proportional to azimuthal vorticity (modulo a phase). The dispersion 
relation follows directly. Equation (Al) gives a simple relation between density and azimuthal 
velocity perturbations, which are nearly in phase since the growth rate is usually small 
compared to — f g k x AU. 
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